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ABSTRACT 

A numerical method for solving the Euler equations for multiblade rotors has 
been developed and some preliminary results are reported. The numerical scheme 
is a combination of several recent methods and algorithm improvements, adapted 
to the particular requirements of rotor-body interactions. A cylindrical basic grid 
has been used to study conventional multiblade helicopter rotors. Test calculations 
have been made for two- and six-blade rotors in hover and for a two-blade rotor 
in forward flight, under transonic tip conditions but without lift. The results show 
good agreement with experimental data. 

1. INTRODUCTION 

To date, computational work has involved a number of finite-difference or 
finite-volume schemes to simulate helicopter rotor flow by solving the potential, Eu- 
ler or Reynolds-averaged Navier-Stokes equations. For most calculations, the com- 
putational domain covers only one blade arbitrarily and the rotor disk is computed 
only partially in the spanwise and circumferential direction. Therefore some phys- 
ical boundary conditions are difficult to implement, and some artificial boundary 
conditions may be required. This work describes general grid topology which covers 
the entire rotor disk, or the partial rotor disk for the hovering case with physical pe- 
riodic boundary conditions. The topology is also suitable for studying rotor-body 
interaction. Moreover, a new Euler solver for helicopter multiblade rotor flow is 
developed. The algorithm can be extended easily to solve the Reynolds-averaged 
Navier-Stokes equations . However, for a satisfactory viscous solution, a finer grid, 
high-order scheme and good turbulence model would be required. It would become 
too expensive to use in a routine production mode. In this paper, the grid topology, 
numerical algorithm, boundary-condition treatment, and verification are described. 
The computed results for nonlifting multiblade rotors in hover and in forward flight 
are compared with experimental data. The blade-to-blade aerodynamic interaction 
for nonlifting rotors is also studied. 

2. GRID SYSTEM 

Most grid topologies used for helicopter rotor flow until now were similar 
to those used for fixed-wing flow. That may be adequate for a one-blade rotor; 
however, those grid topologies are not appropriate or easily extended for multiblade 
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rotors. Indeed, it would seem quite difficult to generate a one-zone structure grid 
for the entire helicopter. The patched grid has been applied successfully in cascade 
(rotor-stator) flows by Rai. 1 In this paper, a similar cylindrical-type patched grid 
topology is proposed for the multiblade rotor and fuselage combination problem, 
which is shown in Fig. 1. The interpolation at the interface between the moving 
grid and the static grid is simplified by using this topology. A rotating C-type or 
O-type fine inner grid can be embedded to wrap the blade with a viscous zone. In 
the viscous zone, either Navier-Stokes or boundary layer equations would be solved. 

The cylindrical rotating outer grid is generated by solving Poisson equations 
with periodic boundary conditions at blade cuts for a one-blade grid set. The 
periodic boundary conditions are Z\ = Zj max and X\ = Xj max — ^ r . The r is 
the radius of each cylinder, X is the arc length along the circumferential direction 
and n is the number of blades. Therefore, a one-blade grid set can be rotated and 
connected to become a multibiade grid set. There are no discontinuities of grid 
slope between each blade grid. The grid set which is generated for a nonlifting 
four-blade rotor is shown in Fig. 2. In this paper, only the Euler calculation for the 
nonlifting rotors is included. The rotating outer grid provides enough resolution for 
the inviscid flow. 

3. INTEGRAL FORMULATION 

In this investigation, the finite- volume approach is adopted for three reasons: 
(a) it is more natural when considering conservation laws, (b) it can avoid metric 
singularities that are usually associated with finite-difference methods, and (c) it 
can capture free-stream conditions without numerical errors. The general form of 
a conservation law applied to an arbitrary grid cell is 2 


/ 

v(t 2 ) 


UdV 


-/ 

V(ti) 


UdV + 



n • F dS dt - 0 


( 1 ) 


where V(t) is the cell volume, n dS(t) is a vector element of surface area with 
outward normal n, U is a conservative variable per unit volume, and F is the flux 
of U per unit area per unit time. Let u and v(t) be the fluid velocity and cell-surface 
element velocity, respectively. The flux F can then be written as 


F = (u - v)U + G (2) 

where the first term is the convection of U relative to the surface element, and G 
stands for the non-convective part of the flux. Although Eq. (l) is general and rigor- 
ous, the relationships between inertial and blade-oriented coordinate systems must 
be considered for rotorcraft applications. Let ro(£), Vn(t).and fl(<) be the position 
vector of the origin, velocity, and angular velocity of the non-inertial coordinate 
system relative to the inertial frame, respectively. Then v can be written as 


V = v r + v c 


( 3 ) 


where 

v r = v 0 (t) + fl(t) x [r - r 0 (t)] (4) 
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and v c (f) is the surface element velocity relative to the non-inertial frame, such as 
blade flapping, lead-lag motion, etc. In the present calculation, ro(t) is set to be 
zero. Also Vo(<) is zero by setting forward flight velocity at far field. Furthermore, 
velocity v c is set to be zero under the assumption that no solution adaptive grid 
is used. The set of conservative variables U , such as density, momentum per unit 
volume, and total internal energy per unit volume is given by the column vector 

( 5 ) 



One can then define the normal flux component F n = n • F, the normal velocity 
components u n = n • u and v n = n • v = 0(<) • (r x n), and the normal relative 
velocity component u' = n • (u - v) = u n — v n . The term of r x n is not a function 
of time and is fixed at each cell for rigid-body rotation. For inviscid flow, the set of 
variables F n is given by the column vector 





pu' 

F n = 

P 

— 

mu' + pn 


E 


eu' + pu n 


where M, P, and E are the normal flux of mass, momentum, and energy, and p is 
the pressure. The surface normal unit vector n can be obtained by 


n = C(wt)n r 


with the rotation matrix 


C(wf) 


cos(wt) 

sin(ut) 

0 


- sin(uit) 0 
cos(u>t) 0 

0 1 


(7a) 


(7b) 


where n r is the cell surface vector represented with respect to the rotation frame 
and is only calculated once. With the rotation matrix(7b), one need not recalcu- 
late coordinates (x, y, z) of the grid points at each time step. The details of the 
evaluation of n r and other grid quantities can be found in Ref. 2. 

In the present work, the Euler equations are written using the absolute vari- 
ables in the inertial frame (that is Eqs. (5)-(6)) for generality. It may be mentioned, 
however, that the formulation can be easily transformed to the rotational frame by 
using the rotation matrix 


then one can write 



0 0 
C 0 
0 1 


U = CU\ 


Fn = CF; 


(8) 

( 9 ) 
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and a source term C T • CU " can be derived from the time differentiation. U' is the 
vector of so-called absolute variables written in the rotational frame. 3-4 can be 
expressed as 


pu' 

C T (m u' + pn) 
eu' + pu n 


( 10 ) 


For the hovering case with this transformation, U" is not a function of time, and 
therefore an acceleration scheme may be adopted to solve this form of governing 
equation(Eq.(10) plus source terms). However for the forward flight, this benefit 
is lost. Furthermore, some care must be taken in order to reproduce uniform flow 
when using this formula. 2 For these reasons, the inertial formulation was considered 
preferable for this paper. 


4. NUMERICAL ALGORITHM 


Several Euler codes for rotors have been developed. 4-10 Most of them 
use central difference operators with user adjustable dissipation coefficients. The 
present code uses Steger- Warming flux vector splitting 11 in the circumferential, 
(or primary flow) direction and Jameson’s dissipative terms 12 in the other two di- 
rections. The use of flux splitting can avoid the “tuning” of explicit dissipation 
coefficients. Furthermore the linear stability analysis shows that the use of flux 
vector splitting combined with upwind differencing in the main flow direction, and 
central differencing in the other two directions, can result in an unconditionally 
stable two-factored algorithm. 13 

In the present approach, the conservation law is applied to cells defined by 
the primary grid. 2 The van Leer MUSCL(monotone upstream-centered scheme for 
the conservation laws) approach is used to evaluate the conservative variables U~ 
or U + at cell surface. 14-16 They are given by 


^j + l/2 ,k,l ~ Uj,k,l + 4>j {Uj,k,l Uj-i'k,l)/2 
Uj+\/2,k,l ~ Uj+l,k,l ~ j + 2,k,l ~ Uj + l,k,l)/2 


(Ha) 

(lib) 


where first- or second-order approximations are obtained by setting <j> = 0 or <f> 
=1, respectively. The geometric quantities related to Fj + - [ , 2kl are evaluated at 
j + 1/2 ,k,l. The resultant two-factor, flux-splitting algorithm can be summarized 
by 

{V + At£ + ){V + Att~)AU n = RHS (12a) 

£+ = b~A + + 6'C - Aic (12b) 

r = ~ 6'B - D t 'n (12c) 

RHS = -VAt(S~(E + ) n + 6+(E-) n + b c v F n - b c .G n - D e U + (12d) 

LA v 


6?A±AV = Af +t/w 




- A ± 




where the quantities D t , D e are the smoothing terms in the crossflow direction, 13 
j E ± — E ± (U Zf ), A ± = A ± (U Zf ) and A U' and AU + are obtained by 


^ j + l/2,k,l ~ ^Uj,k,l + 2 


(12e) 
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Al/+ + = Al>„ - «t(A£/ )+2 , t ,, - Ai/, +w )/2 


( 12 f) 


In Eq. (12), the superscript n is the number of the time steps. 

5. BOUNDARY CONDITIONS 

For inviscid flow, no flow across the wall boundary is implemented by im- 
posing u n = v n on the solid wall. The normal flux component F n for the wall cell 
surface becomes 



m' 


’ 0 1 

F n = 

p 

E 

— 

? ■§ 

1 


(13) 


The only unknown is pressure, which can be determined by solving the normal mo- 
mentum equation (total enthalpy is not constant). The normal momentum equation 
on the wall can be written as 


-p(m - v,). 


Dn\ 

~Dt 


+ />nj. 


D\ i 
~EH 


= ~Pn 


(14a) 


The velocity components required for the normal momentum equation on the wall 
are determined by 

«i •t i = ( 2 u 2 -u 3 ) ( 14 b) 

Ui -tr, = ( 2 u 2 - U 3 ) -t„ ( 14 c) 


and 

v x • ni = u x • ni ( 14 d) 

and tr, are the vectors joining opposite edge midpoints of cells on the blade 
surface, which are normal to the surface vector n x . The subscript 1 represents the 
blade surface, and 2 or 3 are one or two points away from the blade, respectively. 

For the hovering cases, the periodicity is the physical boundary condition 
and can be implemented at the grid cuts as 


Pi — Pjmazi ^1 — Cjmax 

tlj'max = C(-27T,n)u 1 (15) 

where n is the number of blades and subscripts 1 and jmax refer to £ direction. 
There is no additional interpolation required to obtain flow variables at the grid cut 
because the grid set for each blade is generated to satisfy the periodic condition. 
For the far field, nonreflection boundary conditions are used. 12 For the inboard of 
the blades, no flux is allowed through the inboard cylinder. 

6. NUMERICAL RESULTS AND VERIFICATION 

The first test case is a NACA-0012 infinite aspect ratio wing. A very coarse 
C-H type grid is used. There are 49 points in the chordwise direction, 5 points 
in the spanwise direction, and 10 points in the normal-to-the-blade direction. The 
flow is calculated in both the wing-fixed and inertial coordinate system. The far 
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field conditions were either set to = 0.8, a = 0°, or allowed to remain quiescent 
while moving the body-fitted grid at the same Mach number in the negative x 
direction. The pressure contours around the wing surface show excellent agreement 
between these two approaches(Fig.3), though the Mach contours are different by 
a linear transformation as we should expect(Figs. 4(a)-4(b)). Furthermore, the 
instantaneous streamline pattern shows the doublet characteristic of this wing in 
rectilinear motion(Fig. 5) 

The second case is a nonlifting two-blade rotor in hover, which has been in- 
vestigated experimentally by Caradonna and Tung. 17 In this and succeeding cases, 
the equations are solved in the inertial frame. The blade is an untwisted, untapered 
blade with a constant NACA-0012 profile and an aspect ratio of 6. The flow con- 
dition is Mt — 0.52. The pressure comparison is quite good at various locations as 
shown in Fig. 6. This solution was computed using periodic boundary conditions 
(that is n = 2 in Eq. (15)). 

The third case is a six-blade rotor in hover. The relative Mach number 
contours are shown in Figs. 7(a)-(b) (Color Plate I) for Mt = 0.85 and 0.9. The 
sonic line is represented by the white line in this and succeeding figures. Also, 
the violet color represents supersonic flow with respect to a non-inertial coordinate 
system. The different physical characteristics are observed in the tip region. When 
the tip Mach number is increased, the supersonic region on the blade is delocalized 
and connects with the supersonic far field. This effect was also observed in the 
experiment for the same range of Mach number. 18 

The fourth case is a study of the effects of different numbers of blades. For 
comparison of the nonlifting 2-blade and 6-blade rotor, the grid density is essentially 
the same in both cases. The low pressure peak is found to be diminished slightly, 
and the delocalization phenomenon is also slightly delayed by increasing the number 
of blades as shown in Fig. 8. To understand this blade-to-blade interference effect, it 
is useful to consider incompressible flow over parallel infinite cylinders. The distance 
between each neighbor cylinder is 6. By taking leading terms, the complex potential 
W for the uniform flow over infinite cylinders can be approximated as 


n=N 


w 


fia 


_ va 

a 2 Z ■“ b 2 (z - x n — b ') b 2 (z - x n + 6') 

fl — 1 


n'a 2 


n'a 2 


(b - b') 2 (z -x n -(b- b')') (6 + b’) 2 (z -*»-(& + b')') 


n'a 2 


y!a 2 


(b - b') 2 (z -x n + (b- b'Y) {b + b') 2 {z -x n + {b + b'Y) 


+ 


na 


+ 


Ha 


(26 )*(* -*„-(£)) (26)»(* -*„ + (£)) 


z — X, 


where p = U^a 2 , and additional quantities are defined as follows, 

,2 


6 = T’ < 4± >'>' = <sr 
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In these equations, x n is the location of the center of each cylinder, a is the radius 
of a cylinder, n is the number of cylinders and N is a large number. The pressure 
comparison is shown on Fig. 9. Results are presented for b — 5 radii, 10 radii, and 

oo. When the cylinders get closer, the peak of pressure becomes lower. That is the 
same trend as seen in Fig. 8, even though compressibility, numerical viscosity and 
three-dimensional effects are additional factors. 

The fifth case is a nonlifting forward-flight case, studied experimentally by 
Caradonna et al. 19 The experiment tested a stiff, two-bladed, untwisted, untapered 
blade with a constant NACA-0012 profile. The blade was pressure-instrumented at 
r/R = 0.893. Figure 10 compares computed and experimental results for a blade-tip 
Mach number of 0.8 and an advance ratio n of 0.2. Here, pressure coefficients are 
plotted at a single radial location for a variety of rotor azimuthal angles. The grid 
distributed on each blade half surface contains 61 points in the chordwise direction, 
19 points in the spanwise direction, and 46 points normal to the blade. The total 
number of grid points is about 300,000. 

The computational results agree closely with experimental data. The lo- 
cation of the shock for a variety of azimuthal angles is captured very well by the 
present algorithm. The relative Mach number contours on the rotor disk are shown 
in Fig. 11 (Color Plate II). The blue to red represents the low to high Mach num- 
ber region. The supersonic region is represented by violet outside the white line. 
The results show that the shock exists in the advancing side and disappears in the 
retreating side. The contour plots vary with azimuthal angle, which is different 
from the results for hovering cases shown in Fig. 7. 

7. CONCLUDING REMARKS 

A finite-volume, partially flux-vector split Euler code with a cylindrical grid 
topology for multiblade rotor has been developed. The present results of nonlifting 
multiblade rotors show that this code and grid topology are potentially very useful 
for complex multiblade rotor flow. The computed results show good agreement 
with experimental data in either hover or forward flight. In hovering cases, the 
delocalization phenomenon is observed. Moreover, the effects of different numbers 
of blades for the nonlifting multiblade rotor are also examined. When the number 
of blades is increased, the results indicate that the low-pressure peak is reduced 
slightly and consequently the delocalization is delayed slightly. In forward-flight 
cases, the shock can also be captured quite well at a variety of rotor azimuthal 
angles. The lifting cases will be studied in the near future, along with an extension 
to rotor- body interactions. 
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ROTATING OUTER GRID 



Fig. 1 Cylindrical grid topology for rotor-fuselage interaction. 



Fig. 2 Cylindrical grid for a four-blade rotor. 
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Fig. 3 Comparison of pressure contours for a NACA-0012 airfoil with two different 
approaches, a) Wing-fixed coordinate system with far field condition M = 
0.8 and a = 0.0°; b) Inertial coordinate system with moving grid. 




Fig. 4 Comparison of Mach contours for a NACA-0012 airfoil with two different 
approaches, a) Wing-fixed coordinate system with far field condition = 
0.8 and a = 0.0°; b) Inertial coordinate system with moving grid. 



Fig. 5 Instantaneous streamline pattern of a NACA-0012 airfoil in rectilinear mo- 
tion. 
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Fig. 6 Surface pressure results at various locations for a hovering two-blade rotor, 
Mt = 0.52, aspect ratio = 6, nonlifting, untwisted, untapered NACA-0012 
blade. 


Fig. 7 (Color Plate I) follows page 2-2-13. 
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Fig. 8 Comparison of surface pressure results for hovering rotors with different 
numbers of blades, Mt = 0.9, aspect ratio = 6, nonlifting, untwisted, unta- 
pered NACA-0012 blade, a) r/R = 0.90; b) r/R = 0.95; c) r/R = 1.0. 



Fig. 9 Comparison of surface pressure results for inviscid incompressible flows over 
infinite cylinders. 
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Fig. 10 Surface pressure results at the advancing side for a forward-flight two-blade 
rotor, Mt = 0.8, n = 0.2, aspect ratio = 7, nonlifting, untwisted, untapered 


NACA-0012 blade. 


Fig. 11 (Color Plate II) follows page 2-2-13. 
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COLOR PHOTOGRAPH 



(a) 



Fig. 7 Mach contours for a hovering six-blade rotor, aspect ratio = 6, nonlifting, 
untwisted, untapered NACA-0012 blade, a) M T - 0.85; b) M T = 0.9. 


2-2 Color Plate / 





Fig, 11 Mach contours at various azimuthal angles for a forward-flight two-blade 
rotor, Mr 0.8, p = 0.2, aspect ratio = 7, nonlifting, untwisted, untapered 
NACA-0012 blade, a) ^ - 0°; b)'!' = 30°; c)^ = 60°; d)'!' - 90°; e)^ - 
120°; f) ^ = 150°. 

2-2 Color Plate II 
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